#How Do Observers Assess Resolve?
#July 7, 2018
#BJPS Replication 5.R

#This file contains the replication code needed to replicate the analysis from Supplementary Appendix §1.2 (power simulations)
#A huge thanks to Anton Strezhnev for his invaluable expertise!

#To generate the necessary dataframes, either run BJPS Replication 1.R, or load the following saved object:
library(here) #Avoids needing to setwd()
dat3 <- get(load(file="Resolve Conjoint Cleaned No US 070718.RData"))

#Libraries used in the code
library(mlogit)
library(ggplot2)

### Clustered SE function
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) 
  }

## Make mlogit.data file for the mlogit package
dat.complete <- dat[,c("ID", "Round", "Country","RegimeType","Capabilities","Stakes","actor","NewLeader","milService","maleLeader","prevRole","prevOpponent","prevActLdr","curBehavior","Chosen")]
## Identifiers for individual tasks
dat.complete$IDRound <- paste(dat.complete$ID, dat.complete$Round)

## Make mlogit data for the model
dat.mlogit <- mlogit.data(dat.complete, choice="chosen", shape = "long", varying=3:14, alt.var="Country", chid.var="IDRound")

### Start by fitting a multinomial logit just to get baseline utilities
results.baseline <- mlogit(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=dat.mlogit)
summary(results.baseline)
baseline.coefs <- coef(results.baseline)

### Assume a data-generating process where the latent utilities are the same for *every other variable* except for the one under investigation
### Assume the design was fixed, respondent's choices are random
### Then, under a conditional logit choice model, calculate utilities/choice probabilities for each profile in the *actual* set of profiles assigned

## Pull out the model.matrix (corresponding to the coef parameters)
conjoint.design <- model.matrix(results.baseline)

##### Warning: this script takes a while...

##############################
#### Power simulation: capabilities
###############################

#### General parameters + placeholders

capabilities_parameters <- seq(from=0, to=0.2, by = 0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(capabilities_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(capabilities_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of capabilities_parameters
for (param in 1:length(capabilities_parameters)){
  ### The parameter value in the conditional logit model
  gen_param_val <- capabilities_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["CapabilitiesHigh"] <- gen_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A has high capabilities, all else equal
  conjoint.design.capHA <- conjoint.design
  conjoint.design.capHA[,"CapabilitiesHigh"][dat.mlogit$Country == "A"] <- 1
  
  utilities.capHA <- as.vector(exp(conjoint.design.capHA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.capHA <- as.vector(tapply(utilities.capHA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.capHA <- utilities.capHA/sum.utilities.capHA
  
  ### What if state A has low capabilities, all else equal
  conjoint.design.capLA <- conjoint.design
  conjoint.design.capLA[,"CapabilitiesHigh"][dat.mlogit$Country == "A"] <- 0
  
  utilities.capLA <- as.vector(exp(conjoint.design.capLA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.capLA <- as.vector(tapply(utilities.capLA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.capLA <- utilities.capLA/sum.utilities.capLA
  
  amce.A <- choice.prob.capHA[dat.mlogit$Country == "A"] - choice.prob.capLA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B has high capabilities, all else equal
  conjoint.design.capHB <- conjoint.design
  conjoint.design.capHB[,"CapabilitiesHigh"][dat.mlogit$Country == "B"] <- 1
  
  utilities.capHB <- as.vector(exp(conjoint.design.capHB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.capHB <- as.vector(tapply(utilities.capHB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.capHB <- utilities.capHB/sum.utilities.capHB
  
  ### What if leader B has low capabilities, all else equal
  conjoint.design.capLB <- conjoint.design
  conjoint.design.capLB[,"CapabilitiesHigh"][dat.mlogit$Country == "B"] <- 0
  
  utilities.capLB <- as.vector(exp(conjoint.design.capLB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.capLB <- as.vector(tapply(utilities.capLB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.capLB <- utilities.capLB/sum.utilities.capLB
  
  amce.B <- choice.prob.capHB[dat.mlogit$Country == "B"] - choice.prob.capLB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["CapabilitiesHigh",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


cap_sim <- data.frame(amce_true, estimated_power)

##############################
#### Power simulation: stakes
###############################

#### General parameters + placeholders

stakes_parameters <- seq(from=0, to=0.2, by = 0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(stakes_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(stakes_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of stakes_parameters
for (param in 1:length(stakes_parameters)){
  ### The parameter value in the conditional logit model
  gen_param_val <- stakes_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["StakesHigh"] <- gen_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A has high stakes, all else equal
  conjoint.design.stakesHA <- conjoint.design
  conjoint.design.stakesHA[,"StakesHigh"][dat.mlogit$Country == "A"] <- 1
  
  utilities.stakesHA <- as.vector(exp(conjoint.design.stakesHA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.stakesHA <- as.vector(tapply(utilities.stakesHA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.stakesHA <- utilities.stakesHA/sum.utilities.stakesHA
  
  ### What if state A has low stakes, all else equal
  conjoint.design.stakesLA <- conjoint.design
  conjoint.design.stakesLA[,"StakesHigh"][dat.mlogit$Country == "A"] <- 0
  
  utilities.stakesLA <- as.vector(exp(conjoint.design.stakesLA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.stakesLA <- as.vector(tapply(utilities.stakesLA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.stakesLA <- utilities.stakesLA/sum.utilities.stakesLA
  
  amce.A <- choice.prob.stakesHA[dat.mlogit$Country == "A"] - choice.prob.stakesLA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B has high stakes, all else equal
  conjoint.design.stakesHB <- conjoint.design
  conjoint.design.stakesHB[,"StakesHigh"][dat.mlogit$Country == "B"] <- 1
  
  utilities.stakesHB <- as.vector(exp(conjoint.design.stakesHB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.stakesHB <- as.vector(tapply(utilities.stakesHB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.stakesHB <- utilities.stakesHB/sum.utilities.stakesHB
  
  ### What if leader B has low stakes, all else equal
  conjoint.design.stakesLB <- conjoint.design
  conjoint.design.stakesLB[,"StakesHigh"][dat.mlogit$Country == "B"] <- 0
  
  utilities.stakesLB <- as.vector(exp(conjoint.design.stakesLB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.stakesLB <- as.vector(tapply(utilities.stakesLB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.stakesLB <- utilities.stakesLB/sum.utilities.stakesLB
  
  amce.B <- choice.prob.stakesHB[dat.mlogit$Country == "B"] - choice.prob.stakesLB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["StakesHigh",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


stakes_sim <- data.frame(amce_true, estimated_power)

##############################
#### Power simulation: democracy
###############################

#### General parameters + placeholders

democracy_parameters <- seq(from=0, to=-.2, by = -.025) ### Vary the logit parameter from 0 to -.2 (you're getting ~99% power at an effect of .04 which corresponds to logit parameters -.2)
amce_true <- rep(NA, length(democracy_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(democracy_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of democracy_parameters
for (param in 1:length(democracy_parameters)){
  ### The parameter value in the conditional logit model
  dem_param_val <- democracy_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["RegimeTypeDemocracy"] <- dem_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A was a democracy, all else equal
  conjoint.design.democracyA <- conjoint.design
  conjoint.design.democracyA[,"RegimeTypeDemocracy"][dat.mlogit$Country == "A"] <- 1
  conjoint.design.democracyA[,"RegimeTypeMixed"][dat.mlogit$Country == "A"] <- 0 
  
  utilities.democracyA <- as.vector(exp(conjoint.design.democracyA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.democracyA <- as.vector(tapply(utilities.democracyA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.democracyA <- utilities.democracyA/sum.utilities.democracyA
  
  ### What if state A was an autocracy, all else equal
  conjoint.design.autocracyA <- conjoint.design
  conjoint.design.autocracyA[,"RegimeTypeDemocracy"][dat.mlogit$Country == "A"] <- 0
  conjoint.design.autocracyA[,"RegimeTypeMixed"][dat.mlogit$Country == "A"] <- 0 
  
  utilities.autocracyA <- as.vector(exp(conjoint.design.autocracyA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.autocracyA <- as.vector(tapply(utilities.autocracyA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.autocracyA <- utilities.autocracyA/sum.utilities.autocracyA
  
  amce.A <- choice.prob.democracyA[dat.mlogit$Country == "A"] - choice.prob.autocracyA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B was a democracy, all else equal
  conjoint.design.democracyB <- conjoint.design
  conjoint.design.democracyB[,"RegimeTypeDemocracy"][dat.mlogit$Country == "B"] <- 1
  conjoint.design.democracyB[,"RegimeTypeMixed"][dat.mlogit$Country == "B"] <- 0 
  
  utilities.democracyB <- as.vector(exp(conjoint.design.democracyB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.democracyB <- as.vector(tapply(utilities.democracyB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.democracyB <- utilities.democracyB/sum.utilities.democracyB
  
  ### What if state B was an autocracy, all else equal
  conjoint.design.autocracyB <- conjoint.design
  conjoint.design.autocracyB[,"RegimeTypeDemocracy"][dat.mlogit$Country == "B"] <- 0
  conjoint.design.autocracyB[,"RegimeTypeMixed"][dat.mlogit$Country == "B"] <- 0 
  
  utilities.autocracyB <- as.vector(exp(conjoint.design.autocracyB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.autocracyB <- as.vector(tapply(utilities.autocracyB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.autocracyB <- utilities.autocracyB/sum.utilities.autocracyB
  
  amce.B <- choice.prob.democracyB[dat.mlogit$Country == "B"] - choice.prob.autocracyB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    ### Doing clustered standard errors here in line with your estimation, but theoretically, I'm not sure they're necessary 
    ### Given the data-generating process for assignment, the new Abadie, Athey, Imbens, Wooldridge paper I think would suggest
    ### No clustering: https://arxiv.org/abs/1710.02926 - But here it's negligible, and I guess you can go with whatever is more conservative?
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["RegimeTypeDemocracy",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}

dem_sim <- data.frame(amce_true, estimated_power)

##############################
#### Power simulation: actor
###############################

#### General parameters + placeholders

actor_parameters <- seq(from=0, to=-0.2, by = -0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(actor_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(actor_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of actor_parameters
for (param in 1:length(actor_parameters)){
  ### The parameter value in the conditional logit model
  gen_param_val <- actor_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["actorAdversary"] <- gen_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A is an adversary, all else equal
  conjoint.design.adverA <- conjoint.design
  conjoint.design.adverA[,"actorAdversary"][dat.mlogit$Country == "A"] <- 1
  
  utilities.adverA <- as.vector(exp(conjoint.design.adverA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.adverA <- as.vector(tapply(utilities.adverA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.adverA <- utilities.adverA/sum.utilities.adverA
  
  ### What if state A is an ally, all else equal
  conjoint.design.allyA <- conjoint.design
  conjoint.design.allyA[,"actorAdversary"][dat.mlogit$Country == "A"] <- 0
  
  utilities.allyA <- as.vector(exp(conjoint.design.allyA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.allyA <- as.vector(tapply(utilities.allyA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.allyA <- utilities.allyA/sum.utilities.allyA
  
  amce.A <- choice.prob.adverA[dat.mlogit$Country == "A"] - choice.prob.allyA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B is an adversary, all else equal
  conjoint.design.adverB <- conjoint.design
  conjoint.design.adverB[,"actorAdversary"][dat.mlogit$Country == "B"] <- 1
  
  utilities.adverB <- as.vector(exp(conjoint.design.adverB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.adverB <- as.vector(tapply(utilities.adverB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.adverB <- utilities.adverB/sum.utilities.adverB
  
  ### What if leader B is an ally, all else equal
  conjoint.design.allyB <- conjoint.design
  conjoint.design.allyB[,"actorAdversary"][dat.mlogit$Country == "B"] <- 0
  
  utilities.allyB <- as.vector(exp(conjoint.design.allyB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.allyB <- as.vector(tapply(utilities.allyB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.allyB <- utilities.allyB/sum.utilities.allyB
  
  amce.B <- choice.prob.adverB[dat.mlogit$Country == "B"] - choice.prob.allyB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["actorAdversary",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


actor_sim <- data.frame(amce_true, estimated_power)


##############################
#### Power simulation: new leader
###############################

#### General parameters + placeholders

newleader_parameters <- seq(from=0, to=-0.2, by = -0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(newleader_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(newleader_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of newleader_parameters
for (param in 1:length(newleader_parameters)){
  ### The parameter value in the conditional logit model
  gen_param_val <- newleader_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["NewLeaderNew"] <- gen_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if leader A is a new leader, all else equal
  conjoint.design.newleadA <- conjoint.design
  conjoint.design.newleadA[,"NewLeaderNew"][dat.mlogit$Country == "A"] <- 1
  
  utilities.newleadA <- as.vector(exp(conjoint.design.newleadA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.newleadA <- as.vector(tapply(utilities.newleadA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.newleadA <- utilities.newleadA/sum.utilities.newleadA
  
  ### What if leader A is an old leader, all else equal
  conjoint.design.oldleadA <- conjoint.design
  conjoint.design.oldleadA[,"NewLeaderNew"][dat.mlogit$Country == "A"] <- 0
  
  utilities.oldleadA <- as.vector(exp(conjoint.design.oldleadA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.oldleadA <- as.vector(tapply(utilities.oldleadA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.oldleadA <- utilities.oldleadA/sum.utilities.oldleadA
  
  amce.A <- choice.prob.newleadA[dat.mlogit$Country == "A"] - choice.prob.oldleadA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if leader B is a new leader, all else equal
  conjoint.design.newleadB <- conjoint.design
  conjoint.design.newleadB[,"NewLeaderNew"][dat.mlogit$Country == "B"] <- 1
  
  utilities.newleadB <- as.vector(exp(conjoint.design.newleadB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.newleadB <- as.vector(tapply(utilities.newleadB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.newleadB <- utilities.newleadB/sum.utilities.newleadB
  
  ### What if leader B is an old leader, all else equal
  conjoint.design.oldleadB <- conjoint.design
  conjoint.design.oldleadB[,"NewLeaderNew"][dat.mlogit$Country == "B"] <- 0
  
  utilities.oldleadB <- as.vector(exp(conjoint.design.oldleadB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.oldleadB <- as.vector(tapply(utilities.oldleadB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.oldleadB <- utilities.oldleadB/sum.utilities.oldleadB
  
  amce.B <- choice.prob.newleadB[dat.mlogit$Country == "B"] - choice.prob.oldleadB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["NewLeaderNew",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


newlead_sim <- data.frame(amce_true, estimated_power)


##############################
#### Power simulation: military service
###############################

#### General parameters + placeholders

milservice_parameters <- seq(from=0, to=0.2, by = 0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(milservice_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(milservice_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of milservice_parameters
for (param in 1:length(milservice_parameters)){
  ### The parameter value in the conditional logit model
  dem_param_val <- milservice_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["milServiceHigh"] <- dem_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if leader A had long military service, all else equal
  conjoint.design.milHighA <- conjoint.design
  conjoint.design.milHighA[,"milServiceHigh"][dat.mlogit$Country == "A"] <- 1
  conjoint.design.milHighA[,"milServiceSome"][dat.mlogit$Country == "A"] <- 0 
  
  utilities.milHighA <- as.vector(exp(conjoint.design.milHighA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.milHighA <- as.vector(tapply(utilities.milHighA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.milHighA <- utilities.milHighA/sum.utilities.milHighA
  
  ### What if leader A had no military service, all else equal
  conjoint.design.milNoneA <- conjoint.design
  conjoint.design.milNoneA[,"milServiceHigh"][dat.mlogit$Country == "A"] <- 0
  conjoint.design.milNoneA[,"milServiceSome"][dat.mlogit$Country == "A"] <- 0 
  
  utilities.milNoneA <- as.vector(exp(conjoint.design.milNoneA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.milNoneA <- as.vector(tapply(utilities.milNoneA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.milNoneA <- utilities.milNoneA/sum.utilities.milNoneA
  
  amce.A <- choice.prob.milHighA[dat.mlogit$Country == "A"] - choice.prob.milNoneA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if leader B had long military service, all else equal
  conjoint.design.milHighB <- conjoint.design
  conjoint.design.milHighB[,"milServiceHigh"][dat.mlogit$Country == "B"] <- 1
  conjoint.design.milHighB[,"milServiceSome"][dat.mlogit$Country == "B"] <- 0 
  
  utilities.milHighB <- as.vector(exp(conjoint.design.milHighB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.milHighB <- as.vector(tapply(utilities.milHighB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.milHighB <- utilities.milHighB/sum.utilities.milHighB
  
  ### What if leader B had no military service, all else equal
  conjoint.design.milNoneB <- conjoint.design
  conjoint.design.milNoneB[,"milServiceHigh"][dat.mlogit$Country == "B"] <- 0
  conjoint.design.milNoneB[,"milServiceSome"][dat.mlogit$Country == "B"] <- 0 
  
  utilities.milNoneB <- as.vector(exp(conjoint.design.milNoneB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.milNoneB <- as.vector(tapply(utilities.milNoneB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.milNoneB <- utilities.milNoneB/sum.utilities.milNoneB
  
  amce.B <- choice.prob.milHighB[dat.mlogit$Country == "B"] - choice.prob.milNoneB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["milServiceHigh",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}

milserv_sim <- data.frame(amce_true, estimated_power)

##############################
#### Power simulation: gender
###############################

#### General parameters + placeholders

gender_parameters <- seq(from=0, to=0.2, by = 0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(gender_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(gender_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of gender_parameters
for (param in 1:length(gender_parameters)){
  ### The parameter value in the conditional logit model
  gen_param_val <- gender_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["maleLeaderMale"] <- gen_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if leader A was male, all else equal
  conjoint.design.maleA <- conjoint.design
  conjoint.design.maleA[,"maleLeaderMale"][dat.mlogit$Country == "A"] <- 1
  
  utilities.maleA <- as.vector(exp(conjoint.design.maleA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.maleA <- as.vector(tapply(utilities.maleA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.maleA <- utilities.maleA/sum.utilities.maleA
  
  ### What if leader A was female, all else equal
  conjoint.design.femaleA <- conjoint.design
  conjoint.design.femaleA[,"maleLeaderMale"][dat.mlogit$Country == "A"] <- 0
  
  utilities.femaleA <- as.vector(exp(conjoint.design.femaleA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.femaleA <- as.vector(tapply(utilities.femaleA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.femaleA <- utilities.femaleA/sum.utilities.femaleA
  
  amce.A <- choice.prob.maleA[dat.mlogit$Country == "A"] - choice.prob.femaleA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if leader B was male, all else equal
  conjoint.design.maleB <- conjoint.design
  conjoint.design.maleB[,"maleLeaderMale"][dat.mlogit$Country == "B"] <- 1
  
  utilities.maleB <- as.vector(exp(conjoint.design.maleB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.maleB <- as.vector(tapply(utilities.maleB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.maleB <- utilities.maleB/sum.utilities.maleB
  
  ### What if leader B was female, all else equal
  conjoint.design.femaleB <- conjoint.design
  conjoint.design.femaleB[,"maleLeaderMale"][dat.mlogit$Country == "B"] <- 0
  
  utilities.femaleB <- as.vector(exp(conjoint.design.femaleB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.femaleB <- as.vector(tapply(utilities.femaleB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.femaleB <- utilities.femaleB/sum.utilities.femaleB
  
  amce.B <- choice.prob.maleB[dat.mlogit$Country == "B"] - choice.prob.femaleB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["maleLeaderMale",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


male_sim <- data.frame(amce_true, estimated_power)

##############################
#### Power simulation: past actions
###############################

#### General parameters + placeholders

prevact_parameters <- seq(from=0, to=-.3, by = -.025) ### Vary the logit parameter from 0 to -.3 
amce_true <- rep(NA, length(prevact_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(prevact_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of prevact_parameters
for (param in 1:length(prevact_parameters)){
  ### The parameter value in the conditional logit model
  dem_param_val <- prevact_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["prevActLdrSame leader, backed down"] <- dem_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A had a different leader who backed down, all else equal
  conjoint.design.diffbackA <- conjoint.design
  conjoint.design.diffbackA[,"prevActLdrSame leader, backed down"][dat.mlogit$Country == "A"] <- 1
  conjoint.design.diffbackA[,"prevActLdrSame leader, stood firm"][dat.mlogit$Country == "A"] <- 0 
  conjoint.design.diffbackA[,"prevActLdrDifferent leader, backed down"][dat.mlogit$Country == "A"] <- 0 

  
  utilities.diffbackA <- as.vector(exp(conjoint.design.diffbackA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.diffbackA <- as.vector(tapply(utilities.diffbackA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.diffbackA <- utilities.diffbackA/sum.utilities.diffbackA
  
  ### What if state A had a different leader who stood firm, all else equal
  conjoint.design.diffstoodA <- conjoint.design
  conjoint.design.diffstoodA[,"prevActLdrSame leader, backed down"][dat.mlogit$Country == "A"] <- 0
  conjoint.design.diffstoodA[,"prevActLdrSame leader, stood firm"][dat.mlogit$Country == "A"] <- 0 
  conjoint.design.diffstoodA[,"prevActLdrDifferent leader, backed down"][dat.mlogit$Country == "A"] <- 0 
  
  utilities.diffstoodA <- as.vector(exp(conjoint.design.diffstoodA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.diffstoodA <- as.vector(tapply(utilities.diffstoodA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.diffstoodA <- utilities.diffstoodA/sum.utilities.diffstoodA
  
  amce.A <- choice.prob.diffbackA[dat.mlogit$Country == "A"] - choice.prob.diffstoodA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B had a different leader who backed down, all else equal
  conjoint.design.diffbackB <- conjoint.design
  conjoint.design.diffbackB[,"prevActLdrSame leader, backed down"][dat.mlogit$Country == "B"] <- 1
  conjoint.design.diffbackB[,"prevActLdrSame leader, stood firm"][dat.mlogit$Country == "B"] <- 0 
  conjoint.design.diffbackB[,"prevActLdrDifferent leader, backed down"][dat.mlogit$Country == "B"] <- 0 
  
  utilities.diffbackB <- as.vector(exp(conjoint.design.diffbackB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.diffbackB <- as.vector(tapply(utilities.diffbackB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.diffbackB <- utilities.diffbackB/sum.utilities.diffbackB
  
  ### What if state B had a different leader who stood firm, all else equal
  conjoint.design.diffstoodB <- conjoint.design
  conjoint.design.diffstoodB[,"prevActLdrSame leader, backed down"][dat.mlogit$Country == "B"] <- 0
  conjoint.design.diffstoodB[,"prevActLdrSame leader, stood firm"][dat.mlogit$Country == "B"] <- 0
  conjoint.design.diffbackB[,"prevActLdrDifferent leader, backed down"][dat.mlogit$Country == "B"] <- 0  
  
  utilities.diffstoodB <- as.vector(exp(conjoint.design.diffstoodB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.diffstoodB <- as.vector(tapply(utilities.diffstoodB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.diffstoodB <- utilities.diffstoodB/sum.utilities.diffstoodB
  
  amce.B <- choice.prob.diffbackB[dat.mlogit$Country == "B"] - choice.prob.diffstoodB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ##############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
      
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["prevActLdrSame leader, backed down",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


prevAct_sim <- data.frame(amce_true, estimated_power)


##############################
#### Power simulation: previous role initiator/target
###############################

#### General parameters + placeholders

prevrole_parameters <- seq(from=0, to=-0.2, by = -0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(prevrole_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(prevrole_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of prevrole_parameters
for (param in 1:length(prevrole_parameters)){
  ### The parameter value in the conditional logit model
  gen_param_val <- prevrole_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["prevRoleInitiator"] <- gen_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A was previously the initiator, all else equal
  conjoint.design.prevInitA <- conjoint.design
  conjoint.design.prevInitA[,"prevRoleInitiator"][dat.mlogit$Country == "A"] <- 1
  
  utilities.prevInitA <- as.vector(exp(conjoint.design.prevInitA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevInitA <- as.vector(tapply(utilities.prevInitA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevInitA <- utilities.prevInitA/sum.utilities.prevInitA
  
  ### What if state A was previously the target, all else equal
  conjoint.design.prevTarA <- conjoint.design
  conjoint.design.prevTarA[,"prevRoleInitiator"][dat.mlogit$Country == "A"] <- 0
  
  utilities.prevTarA <- as.vector(exp(conjoint.design.prevTarA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevTarA <- as.vector(tapply(utilities.prevTarA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevTarA <- utilities.prevTarA/sum.utilities.prevTarA
  
  amce.A <- choice.prob.prevInitA[dat.mlogit$Country == "A"] - choice.prob.prevTarA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B was previously the initiator, all else equal
  conjoint.design.prevInitB <- conjoint.design
  conjoint.design.prevInitB[,"prevRoleInitiator"][dat.mlogit$Country == "B"] <- 1
  
  utilities.prevInitB <- as.vector(exp(conjoint.design.prevInitB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevInitB <- as.vector(tapply(utilities.prevInitB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevInitB <- utilities.prevInitB/sum.utilities.prevInitB
  
  ### What if state B was previously the target, all else equal
  conjoint.design.prevTarB <- conjoint.design
  conjoint.design.prevTarB[,"prevRoleInitiator"][dat.mlogit$Country == "B"] <- 0
  
  utilities.prevTarB <- as.vector(exp(conjoint.design.prevTarB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevTarB <- as.vector(tapply(utilities.prevTarB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevTarB <- utilities.prevTarB/sum.utilities.prevTarB
  
  amce.B <- choice.prob.prevInitB[dat.mlogit$Country == "B"] - choice.prob.prevTarB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["prevRoleInitiator",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


prevrole_sim <- data.frame(amce_true, estimated_power)

##############################
#### Power simulation: previous opponent adversary/ally
###############################

#### General parameters + placeholders

prevopp_parameters <- seq(from=0, to=0.2, by = 0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(prevopp_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(prevopp_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of prevopp_parameters
for (param in 1:length(prevopp_parameters)){
  ### The parameter value in the conditional logit model
  gen_param_val <- prevopp_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["prevOpponentAdversary"] <- gen_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A was previously fighting an adversary, all else equal
  conjoint.design.prevAdvA <- conjoint.design
  conjoint.design.prevAdvA[,"prevOpponentAdversary"][dat.mlogit$Country == "A"] <- 1
  
  utilities.prevAdvA <- as.vector(exp(conjoint.design.prevAdvA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevAdvA <- as.vector(tapply(utilities.prevAdvA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevAdvA <- utilities.prevAdvA/sum.utilities.prevAdvA
  
  ### What if state A was previously fighting an ally, all else equal
  conjoint.design.prevAllyA <- conjoint.design
  conjoint.design.prevAllyA[,"prevOpponentAdversary"][dat.mlogit$Country == "A"] <- 0
  
  utilities.prevAllyA <- as.vector(exp(conjoint.design.prevAllyA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevAllyA <- as.vector(tapply(utilities.prevAllyA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevAllyA <- utilities.prevAllyA/sum.utilities.prevAllyA
  
  amce.A <- choice.prob.prevAdvA[dat.mlogit$Country == "A"] - choice.prob.prevAllyA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B was previously fighting an adversary, all else equal
  conjoint.design.prevAdvB <- conjoint.design
  conjoint.design.prevAdvB[,"prevOpponentAdversary"][dat.mlogit$Country == "B"] <- 1
  
  utilities.prevAdvB <- as.vector(exp(conjoint.design.prevAdvB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevAdvB <- as.vector(tapply(utilities.prevAdvB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevAdvB <- utilities.prevAdvB/sum.utilities.prevAdvB
  
  ### What if state B was previously fighting an ally, all else equal
  conjoint.design.prevAllyB <- conjoint.design
  conjoint.design.prevAllyB[,"prevOpponentAdversary"][dat.mlogit$Country == "B"] <- 0
  
  utilities.prevAllyB <- as.vector(exp(conjoint.design.prevAllyB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.prevAllyB <- as.vector(tapply(utilities.prevAllyB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.prevAllyB <- utilities.prevAllyB/sum.utilities.prevAllyB
  
  amce.B <- choice.prob.prevAdvB[dat.mlogit$Country == "B"] - choice.prob.prevAllyB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["prevOpponentAdversary",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}


prevopp_sim <- data.frame(amce_true, estimated_power)

##############################
#### Power simulation: mobilized troops
###############################

#### General parameters + placeholders

signal_parameters <- seq(from=0, to=0.2, by = 0.025) ### Vary the logit parameter from 0 to 0.2 
amce_true <- rep(NA, length(signal_parameters)) ### The AMCE implied by the logit parameter
estimated_power <- rep(NA, length(signal_parameters)) ## Store the estimated power 
nIter <- 1000 # Number of iterations per simulation
dat.complete.true <- dat.complete ## Store the true dataset temporarily
pval <- .05 # The p-value that you want power calculations for

### For each value of signal_parameters
for (param in 1:length(signal_parameters)){
  ### The parameter value in the conditional logit model
  dem_param_val <- signal_parameters[param]
  
  ### Make an artificial parameter vector
  artificial_parameters <- baseline.coefs
  ### Change the democracy coefficient
  artificial_parameters["curBehaviorMobilized troops"] <- dem_param_val
  
  ##########################
  ### What's the true AMCE?
  ##########################
  
  ########## State A
  ### What if state A mobilized troops, all else equal
  conjoint.design.mobilizeA <- conjoint.design
  conjoint.design.mobilizeA[,"curBehaviorMobilized troops"][dat.mlogit$Country == "A"] <- 1
  conjoint.design.mobilizeA[,"curBehaviorPublic threat"][dat.mlogit$Country == "A"] <- 0 
  
  utilities.mobilizeA <- as.vector(exp(conjoint.design.mobilizeA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.mobilizeA <- as.vector(tapply(utilities.mobilizeA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.mobilizeA <- utilities.mobilizeA/sum.utilities.mobilizeA
  
  ### What if state A did nothing, all else equal
  conjoint.design.NothingA <- conjoint.design
  conjoint.design.NothingA[,"curBehaviorMobilized troops"][dat.mlogit$Country == "A"] <- 0
  conjoint.design.NothingA[,"curBehaviorPublic threat"][dat.mlogit$Country == "A"] <- 0 
  
  utilities.NothingA <- as.vector(exp(conjoint.design.NothingA%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.NothingA <- as.vector(tapply(utilities.NothingA, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.NothingA <- utilities.NothingA/sum.utilities.NothingA
  
  amce.A <- choice.prob.mobilizeA[dat.mlogit$Country == "A"] - choice.prob.NothingA[dat.mlogit$Country == "A"]
  
  ############# State B
  ### What if state B mobilized troops, all else equal
  conjoint.design.mobilizeB <- conjoint.design
  conjoint.design.mobilizeB[,"curBehaviorMobilized troops"][dat.mlogit$Country == "B"] <- 1
  conjoint.design.mobilizeB[,"curBehaviorPublic threat"][dat.mlogit$Country == "B"] <- 0 
  
  utilities.mobilizeB <- as.vector(exp(conjoint.design.mobilizeB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.mobilizeB <- as.vector(tapply(utilities.mobilizeB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.mobilizeB <- utilities.mobilizeB/sum.utilities.mobilizeB
  
  ### What if state B did nothing, all else equal
  conjoint.design.NothingB <- conjoint.design
  conjoint.design.NothingB[,"curBehaviorMobilized troops"][dat.mlogit$Country == "B"] <- 0
  conjoint.design.NothingB[,"curBehaviorPublic threat"][dat.mlogit$Country == "B"] <- 0 
  
  utilities.NothingB <- as.vector(exp(conjoint.design.NothingB%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities.NothingB <- as.vector(tapply(utilities.NothingB, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  choice.prob.NothingB <- utilities.NothingB/sum.utilities.NothingB
  
  amce.B <- choice.prob.mobilizeB[dat.mlogit$Country == "B"] - choice.prob.NothingB[dat.mlogit$Country == "B"]
  
  #### What's the implied AMCE?
  amce_true[param] <- mean(c(amce.A,amce.B)) ### Simulated in-sample AMCE using a logit model
  
  #####################
  ### What are the actual choice probabilities for the data?
  conjoint.design <- model.matrix(results.baseline)
  
  ## Make the utilities (null ones for this case)
  utilities <- as.vector(exp(conjoint.design%*%artificial_parameters))
  
  ## Calculate choice probs for each task
  sum.utilities <- as.vector(tapply(utilities, dat.mlogit$IDRound, function(x) sum(x))[dat.mlogit$IDRound])
  
  choice.prob <- utilities/sum.utilities
  
  dat.complete.true$choice.prob <- choice.prob

  rejectVec <- rep(0, nIter) ### Vector of rejections per simulation
  ############################################
  ### Do the power calculations - Simulate responses under the conditional logit model
  for (i in 1:nIter){
    sim.dat <- dat.complete.true ## Store a temporary data matrix for the simulation
    sim.chosen <- tapply(sim.dat$choice.prob, sim.dat$IDRound, function(x) sample(c("A","B"), size=1, prob=x))[sim.dat$IDRound] ### Simulate the outcome data
    sim.dat$chosenSimulated <- ifelse(sim.dat$Country == sim.chosen, 1, 0) ### Which one got chosen
    
    ### Run the estimator
    reg.model.sim <- lm(chosenSimulated ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=sim.dat)
    sim.coef.test <- cl(sim.dat, reg.model.sim, as.character(sim.dat$ID)) 
    ### Store if we reject or not (p < .05)
    rejectVec[i] <- ifelse(sim.coef.test["curBehaviorMobilized troops",][4] <= pval, 1, 0)

  }
  ### Store the share of rejections
  estimated_power[param] <- mean(rejectVec)
}

signal_sim <- data.frame(amce_true, estimated_power)

##### Figure 6: Power simulations

#Merge each set of simulations into a single dataframe
power_df <- as.data.frame(rbind(cap_sim, stakes_sim, dem_sim, actor_sim, newlead_sim, milserv_sim, male_sim, prevrole_sim, prevopp_sim, prevAct_sim[1:9,], signal_sim))#Let's keep things at the same length for ease of comparability
power_df$factor <- rep(c("Capabilities", "Stakes", "Regime type", "Relationship with USA", "New leader", "Military service", "Gender", "Role in previous crisis", "Opponent in previous crisis", "Past action", "Costly signal"), each=9)

power_df$factor <- as.factor(power_df$factor)
power_df$factor <- relevel(power_df$factor, ref="Costly signal")
power_df$factor <- relevel(power_df$factor, ref="Past action")
power_df$factor <- relevel(power_df$factor, ref="Opponent in previous crisis")
power_df$factor <- relevel(power_df$factor, ref="Role in previous crisis")
power_df$factor <- relevel(power_df$factor, ref="Gender")
power_df$factor <- relevel(power_df$factor, ref="Military service")
power_df$factor <- relevel(power_df$factor, ref="New leader")
power_df$factor <- relevel(power_df$factor, ref="Relationship with USA")
power_df$factor <- relevel(power_df$factor, ref="Regime type")
power_df$factor <- relevel(power_df$factor, ref="Stakes")
power_df$factor <- relevel(power_df$factor, ref="Capabilities")

#Plot focuses on: high capabilities, highs takes, democracy, adversary, new leader, long military service, male leader, initiator, against adversary, same leader backed down, mobilized troops

qplot(amce_true, estimated_power, data=power_df) + geom_line() + facet_wrap(~ factor, scales="free_x") + labs(x="True AMCE", y="Estimated power") + theme_bw() + geom_hline(aes(yintercept=0.8, alpha=0.5)) + guides(alpha=FALSE)